PD ROM 1

home *** CD-ROM | disk | FTP | other *** search

/ PD ROM 1 / PD ROM Volume I - Macintosh Software from BMUG (1988).iso / Programming / Programming Tools / FORTRAN Routines / RUNKUT.DOC < prev next >

Wrap

Text File | 1985-11-19 | 28.6 KB | 794 lines | [TEXT/ttxt]

1. RUNKUT - AN INITIAL VALUE ORDINARY DIFFERENTIAL EQUATION SOLVER. 1.1. INTRODUCTION RUNKUT is a subroutine that contains three Runge-Kutta initial value ordinary differential equation (IVP ODE) solvers They are: - A fourth order Fehlberg method. - A fifth order Verner method. - An seventh order Verner method. RUNKUT solves coupled IV ODE's of the form: Y' = F(x,Y) There is no limit to the number of coupled equations in a system to be solved other than the amount of memory available. The differential equations are defined in a user-supplied subroutine--more about that later. The subroutine is called using the following FORTRAN call statement: CALL RUNKUT(XA,Y,XB,NEQN,TOLA,TOLR,HSTART,WORK, & IMETH,IERROR,ICOM,FUNC) and must be called from a main "calling" program. 1.2. VARIABLES PASSED TO RUNKUT. The variables passed are explained below and are listed in the order in which they appear in the call statement. The superscripted numbers refer to notes appearing at the end of this section. XA (Real-Input) XA is the starting point of the interval over which integration of the ODE's is to be performed.(1) Y (Real array of size NEQN - Input/Output) On entry, Y must contain the initial values of the Y-variables--that is Y must contain the solutions at the point XA. On exit from RUNKUT, Y will contain the solutions at the end of the specified interval--that is at XB.(2) XB (Real-Input) XB specifies the end of the interval over which integration is to be performed--the point at which solutions to the differential equations are desired.(1,2) NEQN (Integer-Input) NEQN specifies the number of coupled differential equations in the system to be solved. TOLA (Real-Input) TOLA is the absolute error tolerance required on the solution.(3) TOLR (Real-Input) TOLR is the relative error tolerance required on the solution.(3) HSTART (Real-Input/Output) On the very first call to RUNKUT in a calling sequence, HSTART must contain a non-zero value as an initial estimate for the step-size.(1) The actual value used is relatively unimportant as the program will optimize step-sizes to keep the tolerances within specified limits and keep the number of function evaluations to a minimum. It is left to the user's discretion to choose a value that will not require excessive readjustment by the program. On exit from RUNKUT, this variable will contain the last step size used by the program. It is recommended that for subsequent calls to RUNKUT, the step size returned by the program be used as the value for HSTART to solve the equations over the immediately following interval. WORK (Real array of minimum size NEQN*17 - Input) This is a work area made available to RUNKUT by the calling program. IMETH (Integer-Input) Indicates to the solver, the order of the method to be used. IMETH can take on the following values: 1. - Fourth order Fehlberg method 2. - fifth order Verner method 3. - seventh order Verner method While solving a given set of differential equations it is possible to change the order of the method over sucessive intervals simply by setting IMETH to the desired value and resetting ICOM(1) to 0 on each subsequent call to RUNKUT in which a change of order is required.(4) IERROR (Integer-Output) IERROR is an error status indicator returned by RUNKUT after each call. It can have one of the following values: 0 - no error. 1 - the sum of TOLA and TOLR is less than 100 times machine epsilon. This can be caused by setting both TOLA and TOLR to 0, or by specifying values that are too small for the program to handle sucessfully. In this case TOLR is set to a default value. 2 - the problem is too stiff for the method to handle, or there is a discontinuity in the function. It is impossible to proceed in either case.(2) 3 - both conditions 1 and 2 above have occured. IERROR is set to 0 on every entry to RUNKUT. ICOM (Integer array of size 4) ICOM is a communications vector. Each array element is used to pass a specific item of information (as defined below) to and from the subroutine. ICOM(1) (Input) ICOM(1) must be set to 0 on the first call to RUNKUT. ICOM(1) is internally set to a non-zero value by the program after the first call. For subsequent calls to RUNKUT with the same set of equations, ICOM(1) must be left at this non-zero value. It is possible to change the order of the method used (see also IMETH) over sucessive intervals. In this case ICOM(1) must be reset to 0 each time the order of the method is changed.(4) ICOM(2) (Input) A flag that indicates to the program whether to perform checks on the specified values of TOLA and TOLR. 0 - no checking of the error tolerances. 1 - program checks if the sum of TOLA and TOLR is at least 100 times machine epsilon. If not, TOLR is set to a default value and the error indicator IERROR is set to 1. If ICOM(2) is set to 0, the program assumes that it will be receiving non-zero values of at least one of TOLA or TOLR. Failure to check this before entry to RUNKUT could result in fatal program errors. Also, as will be shown in an example, decreasing the error tolerances does not always result in improved solution accuracy. ICOM(3) (Input) A flag that indicates to the program whether to use the default values of mimimum and maximum step size set by the subroutine.(5) 0 - default values of the maximum and minimum step size are used. 1 - allows the user to set values for the maximum and minimum permissible step size. These values are set by including the following FORTRAN common block statement in the main "calling" program: COMMON /CONS/HMIN,HMAX ICOM(4) (Output) Status flag that indicates the presence of round-off error in the solution. 0 - No round-off error. 1 - Severe round-off error possible in answer. ICOM(4) is reset to 0 on entry to RUNKUT. FUNC (Function name) FUNC is replaced by the name of the subroutine that is part of the main user routine in which the differential equations are specified (See section on user supplied routines). The subroutine name must be declared external by the following FORTRAN statement: EXTERNAL [SUBROUTINE NAME]. 1.2.1. NOTES 1. XA could be greater than XB. In other words it is possible to use RUNKUT to solve the differential equations in a negative X-direction given initial values at XA. In that case the step-size will also be negative. HSTART could be specified as a negative value, but that is not essential as the program internally adjusts the algebraic sign of the step-size. 2. If a discontinuity or extreme stiffness is encountered, RUNKUT will return in XB, the value of X closest to the point of discontinuity or stiffness, at which the solutions are still within the specified error bounds. The array Y will then contain solutions at this value of X. 3. If the sum of TOLA and TOLR is less than 100 times machine epsilon (for example if neither were defined in the main program), TOLR is set to a default value of 10(6) times machine epsilon. TOLA is not set to any default value. The program uses a combination of TOLA and TOLR to determine the accuracies of the solutions, and specifying TOLR to a certain value almost always gives better results than specifying TOLA to the same value. 4. Changing the order of the method on every subsequent call to RUNKUT with a given system of equations can lead to excessive computaion times since a large number of constants are re-evaluated each time the order is changed. 5. The default values used by the subroutine are 10(-8) for HMIN, and for HMAX as follows: 0.5 - if the fourth order method is used. 1.0 - if the fifth order method is used. 2.0 - if the seventh order method is used. However, if any of these values are larger than the absolute value of the interval width, HMAX is set to the interval width. It is sometimes necessary to specify HMAX to a small value to keep the step sizes to relatively small values that keep the runge-kutta methods within regions of stability. The subroutine detects the onset of problem "stiffness" or the occurence of a discontinuity in the function by checking if the step size is smaller than HMIN. If the equations are being solved in the negative X-direction, HMIN and HMAX are the smallest and largest permissible step sizes in the negative X-direction respectively. In other words, HMIN and HMAX always define absolute constraints on the step size. 6. Machine epsilon on the IBM PC is very much a function of the compiler used (even with an 8087 math coprocessor installed). For example Microsoft's Fortran-77 compiler (with 8087 support) generates an epsilon of about 10(-16), whereas IBM/Ryan-McFarland's Professional Fortran-77 generates an epsilon of about 10(-20) 1.3. THE USER-SUPPLIED ROUTINES Two routines must be supplied: 1. A main "calling" program that calls RUNKUT. It must also define the initial values of Y at XA, set TOLA and TOLR, set ICOM(1) to 0 on the first call to RUNKUT, define the number of equations in the sytstem to be solved, dimension a WORK array for RUNKUT, and check for the error status as passed back from RUNKUT through IERROR. See the examples at the end of this section. 2. A subroutine containing the system of differential equations. The subroutine name must be declared EXTERNAL in the main calling program and its name passed to RUNKUT through the parameter FUNC-see above. The subroutine is defined using the follwing FORTRAN statement: SUBROUTINE [FUNC] (X,Y,YPRIME,NEQN) where the parameters have the following definitions: X (Real) Contains the current value of X as set in RUNKUT. Do not alter this value within [func] Y (Real - array of size NEQN) Contains the current solutions at X. Do not alter in [func]. YPRIME (Real - array of size NEQN) YPRIME is set in the subroutine to the differentials dy/dx as follows: YPRIME(1) = F/1/(X,Y(1),...,Y(NEQN)) . . YPRIME(NEQN)=F/NEQN/(X,Y(1),...,Y(NEQN)) where F/i/ is the i'th differential function. 1. RUNKUT - AN INITIAL VALUE ORDINARY DIFFERENTIAL EQUATION SOLVER. 1.1. INTRODUCTION RUNKUT is a subroutine that contains three Runge-Kutta initial value ordinary differential equation (IVP ODE) solvers They are: - A fourth order Fehlberg method. - A fifth order Verner method. - An seventh order Verner method. RUNKUT solves coupled IV ODE's of the form: Y' = F(x,Y) There is no limit to the number of coupled equations in a system to be solved other than the amount of memory available. The differential equations are defined in a user-supplied subroutine--more about that later. The subroutine is called using the following FORTRAN call statement: CALL RUNKUT(XA,Y,XB,NEQN,TOLA,TOLR,HSTART,WORK, & IMETH,IERROR,ICOM,FUNC) and must be called from a main "calling" program. 1.2. VARIABLES PASSED TO RUNKUT. The variables passed are explained below and are listed in the order in which they appear in the call statement. The superscripted numbers refer to notes appearing at the end of this section. XA (Real-Input) XA is the starting point of the interval over which integration of the ODE's is to be performed.(1) Y (Real array of size NEQN - Input/Output) On entry, Y must contain the initial values of the Y-variables--that is Y must contain the solutions at the point XA. On exit from RUNKUT, Y will contain the solutions at the end of the specified interval--that is at XB.(2) XB (Real-Input) XB specifies the end of the interval over which integration is to be performed--the point at which solutions to the differential equations are desired.(1,2) NEQN (Integer-Input) NEQN specifies the number of coupled differential equations in the system to be solved. TOLA (Real-Input) TOLA is the absolute error tolerance required on the solution.(3) TOLR (Real-Input) TOLR is the relative error tolerance required on the solution.(3) HSTART (Real-Input/Output) On the very first call to RUNKUT in a calling sequence, HSTART must contain a non-zero value as an initial estimate for the step-size.(1) The actual value used is relatively unimportant as the program will optimize step-sizes to keep the tolerances within specified limits and keep the number of function evaluations to a minimum. It is left to the user's discretion to choose a value that will not require excessive readjustment by the program. On exit from RUNKUT, this variable will contain the last step size used by the program. It is recommended that for subsequent calls to RUNKUT, the step size returned by the program be used as the value for HSTART to solve the equations over the immediately following interval. WORK (Real array of minimum size NEQN*17 - Input) This is a work area made available to RUNKUT by the calling program. IMETH (Integer-Input) Indicates to the solver, the order of the method to be used. IMETH can take on the following values: 1. - Fourth order Fehlberg method 2. - fifth order Verner method 3. - seventh order Verner method While solving a given set of differential equations it is possible to change the order of the method over sucessive intervals simply by setting IMETH to the desired value and resetting ICOM(1) to 0 on each subsequent call to RUNKUT in which a change of order is required.(4) IERROR (Integer-Output) IERROR is an error status indicator returned by RUNKUT after each call. It can have one of the following values: 0 - no error. 1 - the sum of TOLA and TOLR is less than 100 times machine epsilon. This can be caused by setting both TOLA and TOLR to 0, or by specifying values that are too small for the program to handle sucessfully. In this case TOLR is set to a default value. 2 - the problem is too stiff for the method to handle, or there is a discontinuity in the function. It is impossible to proceed in either case.(2) 3 - both conditions 1 and 2 above have occured. IERROR is set to 0 on every entry to RUNKUT. ICOM (Integer array of size 4) ICOM is a communications vector. Each array element is used to pass a specific item of information (as defined below) to and from the subroutine. ICOM(1) (Input) ICOM(1) must be set to 0 on the first call to RUNKUT. ICOM(1) is internally set to a non-zero value by the program after the first call. For subsequent calls to RUNKUT with the same set of equations, ICOM(1) must be left at this non-zero value. It is possible to change the order of the method used (see also IMETH) over sucessive intervals. In this case ICOM(1) must be reset to 0 each time the order of the method is changed.(4) ICOM(2) (Input) A flag that indicates to the program whether to perform checks on the specified values of TOLA and TOLR. 0 - no checking of the error tolerances. 1 - program checks if the sum of TOLA and TOLR is at least 100 times machine epsilon. If not, TOLR is set to a default value and the error indicator IERROR is set to 1. If ICOM(2) is set to 0, the program assumes that it will be receiving non-zero values of at least one of TOLA or TOLR. Failure to check this before entry to RUNKUT could result in fatal program errors. Also, as will be shown in an example, decreasing the error tolerances does not always result in improved solution accuracy. ICOM(3) (Input) A flag that indicates to the program whether to use the default values of mimimum and maximum step size set by the subroutine.(5) 0 - default values of the maximum and minimum step size are used. 1 - allows the user to set values for the maximum and minimum permissible step size. These values are set by including the following FORTRAN common block statement in the main "calling" program: COMMON /CONS/HMIN,HMAX ICOM(4) (Output) Status flag that indicates the presence of round-off error in the solution. 0 - No round-off error. 1 - Severe round-off error possible in answer. ICOM(4) is reset to 0 on entry to RUNKUT. FUNC (Function name) FUNC is replaced by the name of the subroutine that is part of the main user routine in which the differential equations are specified (See section on user supplied routines). The subroutine name must be declared external by the following FORTRAN statement: EXTERNAL [SUBROUTINE NAME]. 1.2.1. NOTES 1. XA could be greater than XB. In other words it is possible to use RUNKUT to solve the differential equations in a negative X-direction given initial values at XA. In that case the step-size will also be negative. HSTART could be specified as a negative value, but that is not essential as the program internally adjusts the algebraic sign of the step-size. 2. If a discontinuity or extreme stiffness is encountered, RUNKUT will return in XB, the value of X closest to the point of discontinuity or stiffness, at which the solutions are still within the specified error bounds. The array Y will then contain solutions at this value of X. 3. If the sum of TOLA and TOLR is less than 100 times machine epsilon (for example if neither were defined in the main program), TOLR is set to a default value of 10(6) times machine epsilon. TOLA is not set to any default value. The program uses a combination of TOLA and TOLR to determine the accuracies of the solutions, and specifying TOLR to a certain value almost always gives better results than specifying TOLA to the same value. 4. Changing the order of the method on every subsequent call to RUNKUT with a given system of equations can lead to excessive computaion times since a large number of constants are re-evaluated each time the order is changed. 5. The default values used by the subroutine are 10(-8) for HMIN, and for HMAX as follows: 0.5 - if the fourth order method is used. 1.0 - if the fifth order method is used. 2.0 - if the seventh order method is used. However, if any of these values are larger than the absolute value of the interval width, HMAX is set to the interval width. It is sometimes necessary to specify HMAX to a small value to keep the step sizes to relatively small values that keep the runge-kutta methods within regions of stability. The subroutine detects the onset of problem "stiffness" or the occurence of a discontinuity in the function by checking if the step size is smaller than HMIN. If the equations are being solved in the negative X-direction, HMIN and HMAX are the smallest and largest permissible step sizes in the negative X-direction respectively. In other words, HMIN and HMAX always define absolute constraints on the step size. 6. Machine epsilon on the IBM PC is very much a function of the compiler used (even with an 8087 math coprocessor installed). For example Microsoft's Fortran-77 compiler (with 8087 support) generates an epsilon of about 10(-16), whereas IBM/Ryan-McFarland's Professional Fortran-77 generates an epsilon of about 10(-20) 1.3. THE USER-SUPPLIED ROUTINES Two routines must be supplied: 1. A main "calling" program that calls RUNKUT. It must also define the initial values of Y at XA, set TOLA and TOLR, set ICOM(1) to 0 on the first call to RUNKUT, define the number of equations in the sytstem to be solved, dimension a WORK array for RUNKUT, and check for the error status as passed back from RUNKUT through IERROR. See the examples at the end of this section. 2. A subroutine containing the system of differential equations. The subroutine name must be declared EXTERNAL in the main calling program and its name passed to RUNKUT through the parameter FUNC-see above. The subroutine is defined using the follwing FORTRAN statement: SUBROUTINE [FUNC] (X,Y,YPRIME,NEQN) where the parameters have the following definitions: X (Real) Contains the current value of X as set in RUNKUT. Do not alter this value within [func] Y (Real - array of size NEQN) Contains the current solutions at X. Do not alter in [func]. YPRIME (Real - array of size NEQN) YPRIME is set in the subroutine to the differentials dy/dx as follows: YPRIME(1) = F/1/(X,Y(1),...,Y(NEQN)) . . YPRIME(NEQN)=F/NEQN/(X,Y(1),...,Y(NEQN)) where F/i/ is the i'th differential function.